Transition Matrix Monte Carlo Reweighting and Dynamics 
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We study an induced dynamics in the space of energy of 
single-spin-fiip Monte Carlo algorithm. The method gives an 
efficient reweighting technique. This dynamics is shown to 
have relaxation times proportional to the specific heat. Thus, 
it is plausible for a logarithmic factor in the correlation time 
of the standard 2D Ising local dynamics. 
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The investigation of new Monte Carlo algorithms has 
been actively pursued due to their importance for effi- 
cient computer simulation. Cluster algorithms [jl]-^ , gen- 
eralized ensemble methods B, together with advanced 
techniques of analysis, such as the histogram methods 
[§-f|, greatly enhanced our ability to do efficient com- 
puter simulations. 

Recently, Swendsen and Li f| proposed a novel algo- 
rithm named transition matrix Monte Carlo (TMMC). 
The method first determines a transition matrix 
W(E\E'), to be defined below, by some sampling tech- 
nique. The invariant state of this matrix gives the 
canonical probability distribution of energy P eq {E) cx 
n(E) exp(— E/ksT), where n(E) is the density of states. 
This method has connections with the broad histogram 
method 0, the projection dynamics E(| , and the canon- 
ical transition probability method ]ll{ |, but its starting 
point and considerations are quite different. 

There are three important aspects in this approach: 
(1) The computation of the matrix can be efficiently im- 
plemented with the iV-fold way (2) the method 
gives a new reweighting technique, efficient and sim- 
ple in comparison with multiple histogram method 
(3) the artificial dynamics is of considerable theoreti- 
cal interest — it could imply a logarithmic correction in 
the two-dimensional Ising model critical relaxation, indi- 
rectly supporting Domany's conjecture Jl|. 

The rest of the paper is organized as follows: we 
first define the transition matrix Monte Carlo (TMMC) 
dynamics. We discuss the reweighting method under 
TMMC. We then describe the steps leading to an exact 
result for W(E\E') for the one-dimensional Ising model. 
We give a differential equation obeyed in the large-size 
limit and indicate its consequence. We give arguments 
for relaxation time in one dimension and relate the relax- 
ation time to specific heat. We present simulation data 
to verify the general results. 



Let us consider a single-spin-flip Glauber dynamics Jji[ ] 
of the Ising model which is described in continuous time 
as 
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where N is the total number of spins, and 
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is the flip rate of site i, which depends on the spin 
value at the site i as well as the values of its nearest 
neighbor spins ov And F t is a flip operator such that 
Fj,P(. . . , <7j, . . .) = P(. . . , — Cj, . . •). The transition ma- 
trix Monte Carlo dynamics is defined, based on the above 
dynamics, by 



d -*§^ = Y,W{E\E')P{E',t), 
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where P(E, t) is the probability of having energy E at 
time t, and 
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Although Eq. (0) and (f§ give the same equilibrium dis- 
tribution of energy, they have totally different dynamics. 

We note that the transition matrix is banded along di- 
agonal. The matrix elements are not independent — the 
column sum is zero due to the conservation of total prob- 
ability, and the row sum satisfies J2e W(E'\E)P eq (E) = 
due to the fact that the equilibrium distribution is 
a stationary distribution. The transition matrix also 
satisfies detailed balance condition, W (E\E')P eq (E r ) = 
W(E'\E)P eq (E), inherited from the detailed balance in 
the original dynamics of spins. The detailed balance con- 
ditions put strong constraint on the matrix elements. For 
example, for any three energies with nonzero transition 
rates among them, we have 

W(E\E")W(E"\E')W(E'\E) = W(E\E')W(E'\E")W(E"\E). 

(5) 
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The matrix elements can be computed by Monte Carlo 
sampling as follows. For a configuration a, we consider 
the number N(a, AE) of cases that energy is changed 
by AE, for the N possible single-spin flips. Then for 
AE ? 0, 



W(E + AE\E) = w(AE)(N(a, AE)) H{a)=E , 



(6) 



where the average is over all configurations having en- 
ergy E, and w(AE) = i(l - tanh(A^/(2fc B T)) for the 
Glauber dynamics. Since the quantity in the angular 
brackets of Eq. @ is independent of temperature (or flip 
rates), once it is determined, we can use it at any tem- 
perature (or with any flip rates) for W(E\E') and conse- 
quently the equilibrium distribution at any temperature. 

The microcanonical averages can be computed in any 
ensembles which have the property that equal energy 
states have equal probability. We used canonical sim- 
ulations at a selection of temperatures, so that the total 
histogram H(E) is approximately flat. Of course, the 
total H(E) obtained by adding results at different T is 
not meaningful. However, it is perfectly correct to add 
statistics to (N(<j, A_E)) £ from equilibrium simulations 
at different temperatures. 

The single histogram method || uses H(E) in a canon- 
ical simulation as an estimate to the equilibrium energy 
distribution P{E). Multiple histogram method attempts 
to overcome the problem of narrow window in the distri- 
bution g. In the present method, adding simulations at 
different T is handled rather naturally. Such a "multi- 
ple histogram" TMMC is extremely simple and effective. 
In connection with the broad histogram method 0, we 
note that when T is set to oo, the detailed balance condi- 
tions reduce to the broad histogram equations J^] . Unfor- 
tunately, the dynamics proposed in Ref. is incorrect 
& 

There are a number of different ways of determin- 
ing the canonical distributions from W(E\E'). We can 
solve the linear equation, J2e' W(E\E')P(E') — 0, or 
use the detailed balance condition W(E\E')P(E') = 
W{E'\E)P(E). The latter seems numerically more 
stable — the equations give us a set of recursion relations 
for P(E). In Fig. |l| we present a calculation of the Ising 
model heat capacity on a 64 by 64 lattice with 25 sim- 
ulations at selected temperatures. The errors are small 
for the whole temperature range. 

The one-dimensional dynamics of spins, Eq. (|l|), can 
be solved analytically [[uj. In particular, the full set of 
eigenvalues of T is known JTJj . We show that the matrix 
W(E\E') can be obtained in closed form in one dimen- 
sion. To begin with, the density of states is n(E) = 2C^ k , 
where k = 0, 1, 2, ... , [L/2\ , assuming periodic bound- 
ary condition. The number of possible states with en- 
ergy E I J = —L + 4fc is twice the number of ways (due 
to overall up-down symmetry) of putting 2k unsatisfied 
bonds among L possible positions. 



The summation over spin states with energy E and 
E' can be expressed in terms of a one-dimensional lat- 
tice gas problem. In one dimension, the flip rate can be 
written as Wi(ai) — (1/2) [l — ^a^ai-i + <7i+i)/2] where 
7 = tanh 2K. Notice that for each final configuration a 
associated with E, contributions to the transition from 
E + = E + 4 J to E come only from each pair of nearest 
neighbors satisfied bonds in configuration a. Thus we 
can write 



n(E+)W(E\E^ 
1 + 7 



= E 



H+l, 



(7) 
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where n.; = and 1 for unsatisfied and satisfied bond, 
respectively. Similar expressions can be written down 
associated with the transition from E~ = E — 4J to 
E. The restricted sums are obtained through a partition 
function in a field as a generating function. We find A^T 
LC^r 2 and A" = 
elements are thus 



LC^_2k~2 an< l the transition matrix 



(fc + l)(2fc+l) 
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(1-7). 
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The diagonal term is computed from the relation 
W k -i,k + W Kk + W k+ i, k = 0, 



(8) 
(9) 

(10) 



and the rest of the elements W k ,w — if \k — k'\ > 1. 

The matrix can not be diagonalized analytically in gen- 
eral. Nevertheless, at zero temperature, 7=1, the eigen- 
values can be obtained explicitly as X k = —2(k+ l)(2k + 
1)/(L — 1), which give us relaxation times as — 1/Afc, with 
the longest relaxation time t = (L — l)/2. 

The TMMC for large system follows a simple and in- 
teresting dynamics. It can be shown rigorously, with the 
method of fi-expansion for master equation [l^ ] , that in 
the large-size scaling limit, the process is described by 
the equation 



dP(x,t') _ d fdP(x,t') 



Of 



dx 



d.r 



+ xP(x,t') 



(11) 



where t' and x are properly scaled time and energy devi- 
ation from equilibrium. 



E - u N 



(12) 



{Ndyi 

and t' — bt with 

b = J im 77X7 E W(E\E)(E - E) 2 , (13) 

JV— »oo Ac N — ' 

E 

where uq is the average energy per spin and c' = ^T 2 c 
is the reduced specific heat per spin. The equation is ob- 
tained by replacing E by x and expanding all the terms 
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in small parameter 1/y/N. Details will be presented else- 
where. 

The continuum limit equation describes a constrained 
random walk. There are two competing effects in the 
current j = —dP/dx — xP; while the first term is the 
usual diffusion, the second term keeps the walker at the 
origin. Equilibrium is obtained when j = 0, giving the 
well-known Gaussian distribution P eq [x) oc e~ x I 2 . 

Equation ( pd| ) can be transformed into a one- 
dimensional quantum harmonic oscillator equation with 
the change of variable P(x,t') = e~ x ^ 4 <p(x,t'). The 
relaxational spectrum is discrete and equally spaced 
(in 1/t). In particular, the relaxation modes are 
e -ni -x / 2 H n (x/V2), where H n is Hermite polynomi- 
als. Translating back to the original time t, the relax- 
ation times are l/(nb) oc c. Applying the general result, 
Eq. ( p"3| ) , to the one-dimensional Ising model, we have 

r n = ^-cosh2K, n = 0,1,2,- ••. (14) 
An 

It is instructive to have an intuitive picture of the 
asymptotic dynamics, which we argue as follows. 

The TMMC is equivalent to the following steps (Algo- 
rithm A). 

1. Do perfect microcanonical simulation, i.e., pick a 
state at random among the degenerate states of 
current energy E. 

2. Do one canonical Monte Carlo move, chosen a site 
at random. 

Repeat step (1) and (2) N times, where N is the number 
of spins in the system. This is one TMMC step (sweep). 

Consider the dynamics at very low temperatures. 
Then only two energies Eq (ground state) and E\ = 
Eq + 4J (first excited states) dominate. Consider K 
such that correlation length (£ oc exp(2K)) is about 
the size of the system. P(E ) oc exp(—E /kBT), 
P{Ex) oc L 2 ex.p((-E Q -U)/k B T), and P (Ex)/ P(E Q ) « 
L 2 exp(-4iT) ps 1. 

We now consider the time scales (from the point of 
view of TMMC) that a transition is made Eq —> E\, and 

Ei — > Eq. 

Let's assume that the system is in its ground state 
Eq. Step (1) actually does nothing; in step (2), a 
kink pair (unsatisfied bonds) is excited with probability 
exp(— 4K). Thus, we need exp(4if) moves to produce a 
kink pair. So the time scale in units of TMMC steps is 
t oc exp(AK) / L oc L. We have used the fact that correla- 
tion length in one dimension is proportional to exp(2K) 
and that correlation length is comparable to system size 
L. The reverse process has a similar scale. This argu- 
ment is consistent with the exact calculation. 

A more general argument fLq ] relating the relaxation 
time to the specific heat can be given, as follows: the 
variance of the energy distribution P(E) is related to 



the specific heat as 5E 2 = cNUbT 2 . Transition matrix 
Monte Carlo moves are random walks in energy confined 
in a region of 5E. Thus the typical time for energy vary- 
ing over 6E is proportional to the square of distance in 
energy. Therefore, time is proportional to 5E 2 , when we 
measure in steps of moves. To get r in sweeps, we must 
divide by N. Let a be the typical time for one step of 
walk in energy, the equation should be 

r oc a(5E 2 /N) = ac k B T 2 = a c', (15) 

We can give a a precise meaning as l/(c'b), where b is 
given in Eq. (|l3|). In one dimension, the unit time a 
diverges. However, this does not happen in dimensions 
d > 1 as T c > 0, and a will be some finite value asymp- 
totically independent of N. 

We check the above results in two dimensions by exten- 
sive Monte Carlo simulation using two different methods: 

(1) the microcanonical/canonical algorithm A described 
above, by computing the time correlation functions; and 

(2) by direct computation of W(E\E') with a Wolff clus- 
ter algorithm ||. The eigenvalues are computed with 
standard package ]19[| . While the first method is re- 
stricted to very small sizes, the second method can apply 
to large systems. In Fig. ||we plot only the inverse eigen- 
values. Very good confirmation of the logL behavior for 
the relaxation time is observed. 

This result has serious implication. In the standard 
single-spin-flip dynamics, if we add between canonical 
Monte Carlo moves microcanonical moves, the resulting 
dynamics is TMMC and still has a residue slowing down 
r ~ log!/ at T c . It is thus difficult to imagine how this 
logarithmic factor can be canceled exactly in the origi- 
nal dynamics. In fact, such logarithmic factor is conjec- 
tured Jl3| , while many numerical computations p0[ did 
not seem to find it explicitly. We think that the Domany 
conjecture is still an open question. 

One of the prediction of Eq. ([ll]) is that the eigenvalues 
of W(E\E') in the large-size limit are equally spaced, this 
is indeed the case for large systems. The eigenf unctions 
associated with these eigen modes are compared with nu- 
merical results. In Fig. [| we plot the analytic results 
(curves) together with numerical values (symbols) from 
exact diagonalization of matrix W(E\E') for a three- 
dimensional Ising model of 16 3 at fegT/J = 6.0. There 
are no adjustable parameters in the comparison except 
an overall normalization. 

In conclusion, the transition matrix Monte Carlo shows 
a novel dynamical behavior. We find an unusual critical 
slowing down in one dimension. For two dimensions and 
higher, we conclude that correlation time is proportional 
to the specific heat. This is in contrast with cluster algo- 
rithms where Li and Sokal j2lj showed that the specific 
heat is only a lower bound to the correlation time. While 
the TMMC artificial dynamics is of theoretical interest, 
the reweighting technique is very useful in practice. 
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FIG. 1. The specific heat of the two-dimensional Ising 
model on a 64 2 lattice by the transition matrix Monte Carlo 
reweighting method. The insert shows the relative error with 
respect to the exact result (obtained numerically based on 
[E2j). The simulations are done at 25 temperatures, each with 
lCr Monte Carlo steps with a single-spin-flip dynamics. 
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FIG. 2. Relaxation time for two-dimensional Ising model at 
T c of the transition matrix Monte Carlo dynamics computed 
from the inverse of the greatest non-zero eigenvalues of matrix 
W(E\E'). The matrices are obtained by a Wolff Monte Carlo 
cluster algorithm with 10 7 to 10 9 cluster flips. 
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FIG. 3. The first four eigenmodes of the transition matrix 
Monte Carlo dynamics. The continuous curves are analytic 
results oc e~ x ^ 2 H n (x/\/2). Symbols are from exact diago- 
nalizations of matrix W(E\E') for a three-dimensional Ising 
model on a 16 3 system at fcsT/J = 6.0. For clarity, only 
every fifth data points are plotted. 



